#!/usr/bin/perl

use Astro::Coord;
use POSIX;
use strict;

# arguments are RA, Dec (of source; sexagesimal), LST (of start time; sexagesimal),
# direction (decimal degrees, 0 is positive el),
# spacing (of pointings, in decimal degrees), number (of pointings), cycles (per pointing),
# name (of pointings)

my ($source_ra,$source_dec,$start_lst,
    $direction_angle,$spacing,$num_pointings,
    $num_cycles,$point_name)=@ARGV;

# we need the ra and dec in turns
my @ra_els=split(/\:/,$source_ra);
my @dec_els=split(/\:/,$source_dec);
my $ra_decimal=($ra_els[0]+$ra_els[1]/60.0+$ra_els[2]/3600.0)*15.0;
my $ra_turns=$ra_decimal/360.0;
my $dec_sign=1;
if ($source_dec=~/^\-/){
    $dec_sign=-1;
}
my $dec_decimal=$dec_sign*$dec_els[0]+$dec_els[1]/60.0+$dec_els[2]/3600.0;
$dec_decimal*=$dec_sign;
my $dec_turns=$dec_decimal/360.0;

# and the LST as well
my @lst_els=split(/\:/,$start_lst);
my $lst_decimal=($lst_els[0]+$lst_els[1]/60.0+$lst_els[2]/3600.0)*15.0;
my $lst_turns=$lst_decimal/360.0;

# the hour angle is just the LST-RA
my $ha_turns=$lst_turns-$ra_turns;
print "starting ha = ".$ha_turns."\n";

# we assume the latitude of the ATCA
my $atca_latitude=-1*(30.0+18.0/60.0+46.385/3600.0);
my $latitude_turns=$atca_latitude/360.0;

# now we get the start az & el
my ($start_az,$start_el)=eqazel($ha_turns,$dec_turns,$latitude_turns,0);

print "starting az,el = ".($start_az*360.0).",".($start_el*360.0)."\n";

open(MOSAIC,">$point_name.mos");
my $ra_offset=0.0;
my $dec_offset=0.0;
my $az_offset=0.0;
my $el_offset=0.0;
my $d_to_r=3.141592654/180.0;
my $cycle_time=10.0*15.0/3600.0; # assume 10 seconds cycle time
for (my $i=0;$i<$num_pointings;$i++){
#    print MOSAIC $ra_offset." ".$dec_offset." ".$num_cycles." \$".$point_name."_".($i+1)."\n";
    print MOSAIC $az_offset." ".$el_offset." -".$num_cycles." \$".$point_name."_".($i+1)."\n";
    # now add the total cycle time to the LST
    $lst_decimal+=$num_cycles*$cycle_time;
    my ($new_lst,$new_lst_hour,$new_lst_minute,$new_lst_second);
    $new_lst=$lst_decimal/15.0;
    $new_lst_hour=floor($new_lst);
    $new_lst-=$new_lst_hour;
    $new_lst*=60.0;
    $new_lst_minute=floor($new_lst);
    $new_lst-=$new_lst_minute;
    $new_lst*=60.0;
    $new_lst_second=floor($new_lst);
    print " new LST = ".$new_lst_hour.":".$new_lst_minute.":".$new_lst_second."\n";
    $lst_turns=$lst_decimal/360.0;
    $ha_turns=$lst_turns-$ra_turns;
    # now figure out what the new az and el have to be
    my $el_inc=$spacing*cos($direction_angle*$d_to_r);
    if (abs($el_inc)<1e-8){
	$el_inc=0;
    }
    $el_offset+=$el_inc;
    my $az_inc=$spacing*sin($direction_angle*$d_to_r);
    if (abs($az_inc)<1e-8){
	$az_inc=0;
    }
    $az_offset+=$az_inc;
    
    print " offsets in az,el = ".$az_offset.",".$el_offset."\n";
    # what will the new reference az,el be?
    my ($ref_az,$ref_el)=eqazel($ha_turns,$dec_turns,$latitude_turns,0);
    print " new ref az,el = ".($ref_az*360.0).",".($ref_el*360.0)."\n";
    # and the appropriate actual az,el we want to look at
    my $act_az=$ref_az+$az_offset/360.0;
    my $act_el=$ref_el+$el_offset/360.0;
    print " new req az,el = ".($act_az*360.0).",".($act_el*360.0)."\n";
    # now get the corresponding ha,dec
    my ($new_ha,$new_dec)=eqazel($act_az,$act_el,$latitude_turns,1);
#    if ($act_az<180.0/360.0){
#	# we're east
#	$new_ha-=1.0;
#    }
    print " new ha = ".($new_ha)."\n";
    # what is that in RA?
    my $new_ra=$lst_turns-$new_ha;
    print " new ra,dec = ".($new_ra*360.0).",".($new_dec*360.0)."\n";
    $ra_offset=($new_ra-$ra_turns)*360.0;
    $dec_offset=($new_dec-$dec_turns)*360.0;
}
close(MOSAIC);
